#!/usr/bin/perl
# 2012/10/17 Get restriction fragment information for each read (filter out >4bp away from restriction site)

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [bowtie aligned file] -o [output file]\n";
}

my %opt;
getopts("i:o:", \%opt);
my $FILE_bowtie = $opt{i};
my $FILE_out = $opt{o};

#---------------------------------------
# read restriction file
#---------------------------------------
my %num2pos;	# fragment numberから位置を示す
my %cate2num;	# categoryからそこに属するfragment numberの一群を示す

#number	chr	position	length_before	length_after
#100001	chr1	1004	1004	230
#100002	chr1	1234	230	1151
#100003	chr1	2385	1151	822
#100004	chr1	3207	822	30
my $FILE_enzyme = '/home/hideki/Project/2012-10-16_DamID_KD_swi6_cut14/out/DpnI_sites.txt';
my $fh_en = IO::File->new($FILE_enzyme) or die "cannot open $FILE_enzyme: $!";
$fh_en->getline();
while($_ = $fh_en->getline()){
	s/\r?\n//;
	my ($number, $chr, $position, $lenB, $lenA) = split /\t/;
	my $cate = int($position / 1000) * 1000;
	$num2pos{$number} = $position;
	push @{$cate2num{$cate}}, $number;
}
$fh_en->close();


#---------------------------------------
# read target file
#---------------------------------------
# 断片の端っこがlocationの位置
#ILLUMINA-CDA390:98:FC:3:1:17901:1020    +       chr3    2099840 1
#ILLUMINA-CDA390:98:FC:3:1:1504:1020     -       chr3    11248   1
#ILLUMINA-CDA390:98:FC:3:1:7858:1020     +       chr3    10890   1
#ILLUMINA-CDA390:98:FC:3:1:15874:1020    -       chr1    4055757 1
#ILLUMINA-CDA390:98:FC:3:1:16453:1021    +       chr3    7572    1
#ILLUMINA-CDA390:98:FC:3:1:6192:1021     -       chr3    2442366 1
#ILLUMINA-CDA390:98:FC:3:1:17211:1021    -       chr2    3012381 1
#ILLUMINA-CDA390:98:FC:3:1:18089:1021    -       chr3    8530    1
#ILLUMINA-CDA390:98:FC:3:1:13771:1021    -       chr3    17556   1


sub searchList{
	my ($cate, $ref_cate, $ref_total) = @_;
	if($$ref_total != 0){
		return;
	}
	if(exists $cate2num{$cate}){
		$ref_cate->{$cate} = 1;
		$$ref_total++;
	}else{
		my $cate_plus = $cate + 1000;
		&searchList($cate_plus, $ref_cate, $ref_total);
		my $cate_minus = $cate - 1000;
		&searchList($cate_minus, $ref_cate, $ref_total);
	}
}

# summary data
my %observed;

my $fh_bowtie = IO::File->new($FILE_bowtie) or die "cannot open $FILE_bowtie: $!";
while($_ = $fh_bowtie->getline()){
	s/\r?\n//;
	my ($id, $direction, $chr, $position, $tmp) = split /\t/;
	my $cate = int($position / 1000) * 1000;

	my %list;
	my $total = 0;
	&searchList($cate, \%list, \$total);

	my @list;
	foreach my $c(keys %list){
		push @list, @{$cate2num{$c}};
	}

	my @sorted = sort {$a <=> $b } @list;
	my $MAX_plus1 = $sorted[$#sorted] + 1;
	my $min_minus1 = $sorted[0] - 1;

	my $min_distance = 9999999;
	my $min_fragment = 0;
	foreach my $f(@sorted, $MAX_plus1, $min_minus1){
		if(exists $num2pos{$f}){
			my $restriction_site = $num2pos{$f};
			my $distance = 0;
			$distance = abs($position - $restriction_site);
			if($distance < $min_distance){
				$min_distance = $distance;
				$min_fragment = $f;
			}
		}
	}

	# もし制限酵素部位から4bp以内の距離だったら情報を取得して出力する
	if($min_distance > 4){
		next;
	}
	my $ID = $min_fragment;
	my $side = '';
	if($direction eq '+'){
		$side = 'R';
	}elsif($direction eq '-'){
		$side = 'L';
	}
	$observed{$ID}{$side}++;
}
$fh_bowtie->close();

#---------------------------------------
# output
#---------------------------------------
# siteNumber chr position L/R(left or right) observed_number
my $fh_out = IO::File->new($FILE_out, 'w') or die "cannot write $FILE_out: $!";
$fh_out->print(join("\t", qw(SiteNumber chr position LR observedNumber)) . "\n");
foreach my $f(sort {$a <=> $b} keys %observed){
	my $position = $num2pos{$f};
	my $chr = '';
	if($f > 300000){
		$chr = 'chr3';
	}elsif($f > 200000){
		$chr = 'chr2';
	}else{
		$chr = 'chr1';
	}
	if(exists $observed{$f}{L}){
		$fh_out->printf("%d\t%s\t%d\t%s\t%d\n", $f, $chr, $position, "L", $observed{$f}{L});
	}
	if(exists $observed{$f}{R}){
		$fh_out->printf("%d\t%s\t%d\t%s\t%d\n", $f, $chr, $position, "R", $observed{$f}{R});
	}
}
$fh_out->close();
